program BuoyRead
!
! Program to retrieve data from a BUFR file containing buoy information
! taken from the ECMWF MARS
! It also uses a file containing buoy meta data and a trusted list of
! buoys, both in ASCII format
! Output data are also written in ASCII format
!
! Last Modified on: $Date: 2016-07-19 11:44:16 +0200 (Tue, 19 Jul 2016) $
!
!
  use BufrMod, only: set_BUFR_fileattributes, get_BUFR_nr_of_messages, &
                     open_BUFR_file, get_BUFR_message, close_BUFR_file, &
                     bufr_file_attr_data, BufrDataType, &
                     max_nr_of_descriptors, missing_indicator_real

  use Compiler_Features, only: getarg_genscat, iargc_genscat

  use windstress, only: calc_parameters_lkb, calc_u_at_new_height_lkb

  implicit none

  real, parameter :: k_to_c = 273.16  ! Kelvin to Celsius conversion

  type :: buoy_metadata_type
    integer :: buoy_id
    real    :: anemometer_height
    real    :: pressure_sensor_height
    real    :: air_temp_sensor_height
    real    :: sea_temp_sensor_height
    logical :: trusted
  end type buoy_metadata_type

  type(buoy_metadata_type), dimension(1000) :: buoy_metadata
  type(bufr_file_attr_data) :: bufr_file_attributes
  type(BufrDataType)        :: tbd
  integer            :: kmsg, ksub, nr_of_messages
  integer            :: kbuoy, nbuoys, kobs_w, nobs_w, subtype
  integer            :: narg, iarg, error_flag
  integer            :: pos_dtg, pos_lat, pos_lon, pos_bid
  integer            :: pos_spd, pos_dir, pos_tos, pos_ahi
  integer            :: pos_tem, pos_sst, pos_prs, pos_rhu
  integer            :: bid
  real               :: lat, lon
  real               :: spd, dir, tos, ahi, tem, sst, prs, rhu
  real               :: u_star, z0, zl, Cm
  logical            :: duplicate
  character(len=256) :: filename_bufr, filename_output
  character(len=256) :: filename_metadata, filename_trusted_list
  character(len=256) :: arg, arg2, BUFR_TABLES, tmpstr
  character(len=20)  :: dtg
  character(len=100) :: line_to_write

  logical            :: neutral_stability = .true.

  ! keep list of written data to check for double entries
  character(len=45), dimension(20000) :: buoy_data_list

  integer,parameter  :: unit_asc = 29

  ! in this program, we convert all data from the BUFR module to default
  ! size reals (4 bytes, generally)
  ! to prevent rounding errors when checking for missing values, the
  ! value of missing_indicator_real is also rounded
  real,parameter     :: missing_real = missing_indicator_real

  call GetEnv('BUFR_TABLES', BUFR_TABLES)   ! Environment variable
  if (len_trim(BUFR_TABLES) .eq. 0) then
    print '(A)','Please set ${BUFR_TABLES}! Processing halted!'
    stop
  endif

  ! read command line
  narg = iargc_genscat()
  if (narg .eq. 0) then
    print '(A)',' '
    print '(A)','Usage: BuoyRead -f <BUFR input file> -o <ASCII output file>'
    print '(A)','                -bm <buoy metadata file> -tl <trusted list file>'
    print '(A)','with additional option -stab to include stability'
    print '(A)','default: neutral stability'
    print '(A)',' '
    stop
  endif

  iarg = 0
  do while (iarg .lt. narg)
    iarg = iarg + 1
    call getarg_genscat(iarg,arg)

    arg2 = arg(2:)
    select case(arg2)

    ! input file
    case('f')
      iarg = iarg + 1
      call getarg_genscat(iarg,filename_bufr)

    ! output file
    case('o')
      iarg = iarg + 1
      call getarg_genscat(iarg,filename_output)

    ! buoy metadata file
    case('bm')
      iarg = iarg + 1
      call getarg_genscat(iarg,filename_metadata)

    ! list containing trusted buoy id's
    case('tl')
      iarg = iarg + 1
      call getarg_genscat(iarg,filename_trusted_list)

    ! optional inclusion of stability
    case ('stab')
      neutral_stability = .false.
    end select
  enddo

  ! read list of buoy metadata
  error_flag = 0
  nbuoys = 0
  open(unit_asc, file=trim(filename_metadata), status='old')
  do while (error_flag .eq. 0 .and. nbuoys .lt. size(buoy_metadata))
    read(unit_asc, fmt='(A)', iostat=error_flag) tmpstr
    if (error_flag .eq. 0) then
      nbuoys = nbuoys + 1
      buoy_metadata(nbuoys)%trusted = .false.
      read(tmpstr, fmt='(I5,4F6.1)', iostat=error_flag) buoy_metadata(nbuoys)%buoy_id, &
        buoy_metadata(nbuoys)%anemometer_height, buoy_metadata(nbuoys)%pressure_sensor_height, &
        buoy_metadata(nbuoys)%air_temp_sensor_height, buoy_metadata(nbuoys)%sea_temp_sensor_height

      ! buoy id should contain only digits, otherwise skip it
      if (error_flag .ne. 0) then
        nbuoys = nbuoys - 1
        error_flag = 0
      endif

      ! last three digits ge 500 indicates a drifter, skip it from the list
      if (mod(buoy_metadata(nbuoys)%buoy_id, 1000) .ge. 500) then
        nbuoys = nbuoys - 1
      endif

      ! moored buoys have anemometer hight between 3 and 10 meters, skip others
      if (buoy_metadata(nbuoys)%anemometer_height .lt. 3.0 .or. &
          buoy_metadata(nbuoys)%anemometer_height .gt. 10.0) then
        nbuoys = nbuoys - 1
      endif
    endif
  enddo
  close(unit_asc)

  ! read list of trusted buoys and set indicator when a buoy is in the list
  error_flag = 0
  open(unit_asc, file=trim(filename_trusted_list), status='old')
  do while (error_flag .eq. 0)
    read(unit_asc, fmt='(A)', iostat=error_flag) tmpstr
    if (error_flag .eq. 0) then
      read(tmpstr, fmt='(I5)', iostat = error_flag) bid

      ! buoy id should contain only digits
      if (error_flag .eq. 0) then
        do kbuoy = 1,nbuoys
          if (bid .eq. buoy_metadata(kbuoy)%buoy_id) then
            buoy_metadata(kbuoy)%trusted = .true.
            exit ! loop over metadata list
          endif
        enddo
      endif
    endif
  enddo
  close(unit_asc)

  ! open BUFR input file and obtain number of messages present
  call set_BUFR_fileattributes(filename_bufr, bufr_file_attributes, error_flag)
  call open_BUFR_file(bufr_file_attributes, error_flag)
  nr_of_messages = get_BUFR_nr_of_messages(bufr_file_attributes)

  ! open ASCII output file and write heading
  open(unit_asc, file=filename_output)
  write(unit_asc,'(A)') &
    '#Buoy    Lat     Lon AnHi TsHi      Date Time AirTem    SST SrfPres RHu Speed    Dir Sp10m'
  nobs_w = 0

  ! loop over messages in BUFR file
  subtype = -1
  msgloop:do kmsg= 1, nr_of_messages
    call get_BUFR_message(bufr_file_attributes, kmsg, tbd, error_flag)

    ! determine type of message depending on BUFR descriptor

    ! -----------------------------------------------------------------------
    ! 308003: RDB data subtype 21
    ! -----------------------------------------------------------------------
    if (tbd%ktdlst(1) .eq. 308003) then

      ! position of relevant BUFR data descriptors
      if (subtype .ne. 21) then
        print '(A)','RDB data subtype 21'
        subtype = 21
      endif
      pos_bid = 1
      pos_tos = 4
      pos_dtg = 5
      pos_lat = 10
      pos_lon = 11
      pos_spd = 17 ! wind speed
      pos_dir = 16 ! wind direction
      pos_ahi = 33 ! anemometer height
      pos_tem = 18 ! air temperature
      pos_sst = 32 ! sea/water temperature
      pos_prs = 13 ! pressure reduced to mean sea level
      pos_rhu = 20 ! relative humidity

    ! -----------------------------------------------------------------------
    ! 308004: RDB data subtype 13
    ! -----------------------------------------------------------------------
    else if (tbd%ktdlst(1) .eq. 308004) then

      ! position of relevant BUFR data descriptors
      if (subtype .ne. 13) then
        print '(A)','RDB data subtype 13'
        subtype = 13
      endif
      pos_bid = 1
      pos_tos = 4
      pos_dtg = 5
      pos_lat = 10
      pos_lon = 11
      pos_spd = 17 ! wind speed
      pos_dir = 16 ! wind direction
      pos_ahi = 34 ! anemometer height
      pos_tem = 18 ! air temperature
      pos_sst = 32 ! sea/water temperature
      pos_prs = 13 ! pressure reduced to mean sea level
      pos_rhu = 20 ! relative humidity

    ! -----------------------------------------------------------------------
    ! unknown format
    ! -----------------------------------------------------------------------
    else
      print '(A,I7)','Unknown buoy BUFR data format, first table D data descriptor =', &
        tbd%ktdlst(1)
      cycle msgloop ! go to next message
    endif

    ! loop over subsets in BUFR message
    do ksub = 1, tbd%nsubsets

      if (tbd%ktdlst(1) .eq. 308003) then

        ! subtype 21, station identifier is numeric
        bid = nint(tbd%values(max_nr_of_descriptors*(ksub-1)+pos_bid))
      else

        ! subtype 13, CCITTIA5 station identifier is contained in cvals character element
        read(tbd%cvals(max_nr_of_descriptors*(ksub-1)+pos_bid)(1:8),'(I8)', iostat = error_flag) bid
        if (error_flag .ne. 0) cycle     ! buoy id contains non-digits indicating this is a ship
      endif

      if (bid .gt. 99999) cycle          ! buoy id should not contain more than 5 digits
      if (mod(bid,1000) .ge. 500) cycle  ! last three digits ge 500 indicates a drifter

      tos = real(tbd%values(max_nr_of_descriptors*(ksub-1)+pos_tos)) ! type of station
      lat = real(tbd%values(max_nr_of_descriptors*(ksub-1)+pos_lat)) ! latitude
      lon = real(tbd%values(max_nr_of_descriptors*(ksub-1)+pos_lon)) ! longitude
      spd = real(tbd%values(max_nr_of_descriptors*(ksub-1)+pos_spd)) ! wind speed
      dir = real(tbd%values(max_nr_of_descriptors*(ksub-1)+pos_dir)) ! wind direction
      ahi = real(tbd%values(max_nr_of_descriptors*(ksub-1)+pos_ahi)) ! anemometer height
      tem = real(tbd%values(max_nr_of_descriptors*(ksub-1)+pos_tem)) ! air temp (K)
      sst = real(tbd%values(max_nr_of_descriptors*(ksub-1)+pos_sst)) ! SST (K)
      prs = real(tbd%values(max_nr_of_descriptors*(ksub-1)+pos_prs)) ! surface pressure
      rhu = real(tbd%values(max_nr_of_descriptors*(ksub-1)+pos_rhu)) ! relative humidity

      if (tos .ne. 0.0) cycle                   ! type of station should be 'automatic'
      if (lat .eq. missing_real) cycle
      if (lon .eq. missing_real) cycle
      if (spd .le. 0.1) cycle                   ! wind speed QC
      if (spd .ge. 50.0) cycle                  ! wind speed QC
      if (dir .eq. missing_real) cycle
      if (tem .le. 233.16) cycle                ! air temp below -40 C
      if (tem .ge. 313.16) cycle                ! air temp above 40 C
      if (sst .le. 271.16) cycle                ! SST below -2 C
      if (sst .ge. 313.16) cycle                ! SST above 40 C
      if (prs .le. 84000.0) prs = -1.0          ! pressure not mandatory, set it to missing
      if (prs .ge. 105000.0) prs = -1.0         ! pressure not mandatory, set it to missing
      if (rhu .le. 0.0) rhu = -1.0              ! rel humidity not mandatory, set it to missing
      if (rhu .ge. 100.0) rhu = -1.0            ! rel humidity not mandatory, set it to missing

      ! find buoy in buoy metadata list
      do kbuoy = 1,nbuoys
        if (bid .eq. buoy_metadata(kbuoy)%buoy_id) then

          ! check anemometer heights in BUFR and buoy metadata list
          ! BUFR anemometer height is sometimes set to missing or 0 but this data is still valid
          if (buoy_metadata(kbuoy)%trusted) then
            if (ahi .lt. missing_real*0.99 .and. ahi .gt. 0.01) then
              if (abs(ahi - buoy_metadata(kbuoy)%anemometer_height) .gt. 0.5) then
                print '(A,I6,2(A,F5.1))', 'Discrepancy in anemometer heights for buoy', bid, &
                  ': in BUFR:', ahi, ', in metadata:', buoy_metadata(kbuoy)%anemometer_height
                buoy_metadata(kbuoy)%trusted = .false.
              endif
            endif
          endif

          ! write data if buoy is trusted
          if (buoy_metadata(kbuoy)%trusted) then

            ! compute LKB stress model parameters, depending on the availability of
            ! relative humidity and surface pressure
            if (rhu .gt. 0.0 .and. prs .gt. 0.0) then
              call calc_parameters_lkb(height_u=buoy_metadata(kbuoy)%anemometer_height, &
                                       height_t=buoy_metadata(kbuoy)%air_temp_sensor_height, &
                                       height_q=buoy_metadata(kbuoy)%air_temp_sensor_height, &
                                       u=spd, sst=sst-k_to_c, t=tem-k_to_c, &
                                       r_in=rhu/100.0, p_in=prs/100.0, &
                                       u_star=u_star, z0=z0, zl=zl, Cm=Cm)
            elseif (rhu .lt. 0.0 .and. prs .gt. 0.0) then
              call calc_parameters_lkb(height_u=buoy_metadata(kbuoy)%anemometer_height, &
                                       height_t=buoy_metadata(kbuoy)%air_temp_sensor_height, &
                                       height_q=buoy_metadata(kbuoy)%air_temp_sensor_height, &
                                       u=spd, sst=sst-k_to_c, t=tem-k_to_c, p_in=prs/100.0, &
                                       u_star=u_star, z0=z0, zl=zl, Cm=Cm)
            elseif (rhu .gt. 0.0 .and. prs .lt. 0.0) then
              call calc_parameters_lkb(height_u=buoy_metadata(kbuoy)%anemometer_height, &
                                       height_t=buoy_metadata(kbuoy)%air_temp_sensor_height, &
                                       height_q=buoy_metadata(kbuoy)%air_temp_sensor_height, &
                                       u=spd, sst=sst-k_to_c, t=tem-k_to_c, r_in=rhu/100.0, &
                                       u_star=u_star, z0=z0, zl=zl, Cm=Cm)
            else
              call calc_parameters_lkb(height_u=buoy_metadata(kbuoy)%anemometer_height, &
                                       height_t=buoy_metadata(kbuoy)%air_temp_sensor_height, &
                                       height_q=buoy_metadata(kbuoy)%air_temp_sensor_height, &
                                       u=spd, sst=sst-k_to_c, t=tem-k_to_c, &
                                       u_star=u_star, z0=z0, zl=zl, Cm=Cm)
            endif

            if (u_star .gt. missing_indicator_real * 0.99) exit
            if (z0     .gt. missing_indicator_real * 0.99) exit
            if (zl     .gt. missing_indicator_real * 0.99) exit
            if (Cm     .gt. missing_indicator_real * 0.99) exit

            ! convert to neutral winds by setting the stability parameter to 0
            if (neutral_stability) zl = 0.0

            write(dtg,101) &   ! write date/time information to variable
                  int(tbd%values(max_nr_of_descriptors*(ksub-1)+pos_dtg)),   &
                  int(tbd%values(max_nr_of_descriptors*(ksub-1)+pos_dtg+1)), &
                  int(tbd%values(max_nr_of_descriptors*(ksub-1)+pos_dtg+2)), &
                  int(tbd%values(max_nr_of_descriptors*(ksub-1)+pos_dtg+3)), &
                  int(tbd%values(max_nr_of_descriptors*(ksub-1)+pos_dtg+4))

            ! write buoy information, measured wind and 10m wind
            write(line_to_write,102) bid, lat, lon, &
              buoy_metadata(kbuoy)%anemometer_height, &
              buoy_metadata(kbuoy)%air_temp_sensor_height, &
              trim(dtg), tem, sst, nint(prs), nint(rhu), spd, dir, &
              calc_u_at_new_height_lkb(buoy_metadata(kbuoy)%anemometer_height, 10.0, spd, z0, zl, Cm)

            ! check if this observation has been written before and prevent duplicates
            duplicate = .false.
            do kobs_w = 1,nobs_w
              if (buoy_data_list(kobs_w) .eq. line_to_write(1:45)) then
                duplicate = .true.
                exit
              endif
            enddo
            if (.not. duplicate) then
              write(unit_asc,'(A)') trim(line_to_write)
              nobs_w = nobs_w + 1
              buoy_data_list(nobs_w) = line_to_write(1:45)
            endif

          endif
          exit ! loop over metadata list
        endif
      enddo

    enddo
  enddo msgloop

  ! close input and output files
  call close_BUFR_file(bufr_file_attributes, error_flag)
  close(unit_asc)

  ! format statements
  101 format (I4.4,2I2.2,1X,2I2.2)
  102 format (I5,F7.2,F8.2,F5.1,F5.1,2X,A,F7.1,F7.1,I8,I4,F6.1,F7.1,F6.1)

end program BuoyRead
